/*
 * genomeindexer.cpp
 *
 *  Created on: Jul 29, 2011
 *      Author: Pierre Lindenbaum PhD
 *              plindenbaum@yahoo.fr
 *              http://plindenbaum.blogspot.com
 *              
 */
#include <string>
#include <stdexcept>
#include <cstdlib>
#include <cstdio>
#include <climits>
#include <cstring>
#include <iostream>
#include <vector>
#include <list>
#include <cerrno>
#include <sys/types.h>
#include <dirent.h>
#include <algorithm>
#include "xzlib.h"
#include "xstdio.h"
#include "xstdlib.h"
#include "xstring.h"
#include "xfaidx.h"
#include "mergesort.h"
#include "indexedpositions.h"
#include "typedfile.h"
#include "bio.h"
#include "where.h"
#include "shortread.h"

//#define TEST_SEQ "GAGACCAACCGAACCCCCTTCGACCTTGCCGAAGGGGAGTCCGAACTAGTCTCAGGCTTCAACATCGAAT"


using namespace std;

typedef int seq_length_t ;



class IndexedPositions:public TypedFile<IndexedPosition>
    {
    private:
	 IndexedGenome* genome;

    public:
	IndexedPositions(IndexedGenome* genome,const char* indexName):
	    TypedFile<IndexedPosition>(indexName),
	    genome(genome)
	    {

	    }


	virtual ~IndexedPositions()
	    {

	    }

	IndexedGenome* getGenome()
	    {
	    return genome;
	    }

	struct AbstractHit
	    {
	    IndexedPosition position;
	    CharSequence const* seq;
	    size_t count;
	    };

	class AlgoFind
	    {
	    protected:

		IndexedPositions* owner;
		virtual AbstractHit* createHit()=0;

		virtual size_t prefixLength(
		    IndexedPosition* pos,
		    const CharSequence* shortRead
		    )
		     {
		     FaidxPtr faidx= owner->genome->getFaidByTid(pos->chrom_id);
		     ChromosomeString chrom(owner->genome,faidx);
		     int64_t position=pos->position;
		     size_t seqidx=0;
		     while(seqidx < shortRead->size() && position< chrom.size())
			  {
			  int c1= normalize_base(chrom.at(position));
			  int c2= normalize_base(shortRead->at(seqidx));
			  int i=c1-c2;
			  if(i!=0) return seqidx;
			  seqidx++;
			  position++;
			  }
		     return seqidx;
		     }

		virtual int compare(IndexedPosition* pos,const CharSequence* seq)
		     {
		     FaidxPtr faidx= owner->genome->getFaidByTid(pos->chrom_id);
		     int64_t position=pos->position;
		     ChromosomeString chrom(owner->genome,faidx);
		     size_t seqidx=0;

		     for(;;)
			  {
			  if(seqidx == seq->size())
				{
				return 0;
				}
			   else if(position >= chrom.size())
				{
				return 1;
				}

			    int c1= normalize_base(chrom.at(position));
			    int c2= normalize_base(seq->at(seqidx));

			    int i=c1-c2;
			    if(i!=0)
				{
				return i;
				}
			    seqidx++;
			    position++;
			    }

		     }

		virtual long long lowerBound(const CharSequence* seq)
		    {
		    IndexedPosition position;
		    long beg=0L;
		    long end=owner->size();
		    long len=end-beg;
		    while(len>0)
			{
			long half = len/2;
			long middle=beg+half;

			owner->get(middle,&position);

			if(compare(&position,seq)<0)
			    {
			    beg = middle;
			    ++beg;
			    len = len - half - 1;
			    }
			 else
			    {
			    len = half;
			    }
			}
		    return beg;
		    }


		virtual bool startsWith(
		    IndexedPosition* pos,
		    const char* seq,size_t seqLen,
		    size_t maxLen
		    )
		     {
		     FaidxPtr faidx= owner->genome->getFaidByTid(pos->chrom_id);
		     int64_t position=pos->position;
		     size_t seqidx=0;
		     for(;;)
			  {
			  if(seqidx == seqLen || seqidx==maxLen)
				{
				return true;
				}
			   else if(position >= faidx->len)
				{
				return false;
				}

			    int c1= normalize_base(owner->genome->at(faidx,position));
			    int c2= normalize_base(seq[seqidx]);
			     int i=c1-c2;
			    if(i!=0) return false;
			    seqidx++;
			    position++;
			    }
		     }
	     public:
		void* userData;
		AlgoFind(IndexedPositions* owner):owner(owner),userData(NULL)
		    {

		    }
		virtual ~AlgoFind()
		    {

		    }

		virtual int64_t find(const CharSequence* seq)=0;
		virtual void hitFound(const AbstractHit* hit){}
	    };

	class PerfectMatch:public AlgoFind
	    {
	    public:
		struct PerfectHit:public AbstractHit
		    {

		    };
	    protected:

		virtual AbstractHit* createHit()
		    {
		    return new PerfectHit;
		    }
	    public:
		PerfectMatch(IndexedPositions* owner):AlgoFind(owner)
		    {

		    }
		virtual ~PerfectMatch() {}

		virtual int compare2(IndexedPosition* pos,const CharSequence* seq,AbstractHit* hit)
		     {
		     FaidxPtr faidx= owner->genome->getFaidByTid(pos->chrom_id);
		     int64_t position=pos->position;
		     size_t seqidx=0;
		     //bool debug=strncasecmp( TEST_SEQ, seq,seqLen)==0;
		     for(;;)
			  {
			  if(seqidx == seq->size())
				{
				return 0;
				}
			   else if(position >= faidx->len)
				{
				return 1;
				}

			    int c1= normalize_base(owner->genome->at(faidx,position));
			    int c2= normalize_base(seq->at(seqidx));
			    //if(debug) cerr << "comparing["<< seqidx<<"/" << seqLen <<"|"<< position <<"]"<<(char)c1 <<" / " << (char)c2 << endl;

			    int i=c1-c2;
			    if(i!=0)
				{
				//if(debug) cerr << "break\n";
				return i;
				}
			    seqidx++;
			    position++;
			    }
		     }

		virtual long long __deprecated_lowerBound(const CharSequence* seq,AbstractHit* hit)
		    {
		    IndexedPosition position;
		    long beg=0L;
		    long end=owner->size();
		    long len=end-beg;
		    //bool debug=strncasecmp( TEST_SEQ, seq,seqLen)==0;
		    //if(debug) cerr << "LEN="<< len << endl;
		    while(len>0)
			{
			long half = len/2;
			long middle=beg+half;

			owner->get(middle,&position);
			//if(debug) cerr << beg << "<" <<  middle << "<" << end << " pos:"<< position.position << endl;
			if(compare2(&position,seq,hit)<0)
			    {
			    beg = middle;
			    ++beg;
			    len = len - half - 1;
			    }
			 else
			    {
			     len = half;
			    }
			}
		    //cerr << "exit loop" << "\n";
		    return beg;
		    }

		virtual int64_t find(const CharSequence* seq)
		    {
		    AbstractHit* hit=createHit();
		    hit->seq=seq;
		    hit->count=0;
		    //bool debug=strncasecmp( TEST_SEQ, seq,lenSeq)==0;
		    //cerr << "OK" << debug<< endl;
		    //cerr << TEST_SEQ << "\n";
		    //cerr.write(seq,lenSeq);
		    //cerr << endl;
		    long n= this->lowerBound(seq);
		    //if(debug) cerr << n << endl;
		    while(n<  owner->size()  )
			{

			owner->get(n,&(hit->position));

			if(this->compare2(&(hit->position),seq,hit)!=0) break;
			hit->count++;
			hitFound(hit);
			++n;
			}
		    int64_t count= hit->count;
		    delete hit;
		    return count;
		    }

	    };



	    class MisMatch:public PerfectMatch
		    {
		    public:
			struct MismatchHit:public AbstractHit
			    {

			    };
			 int maxMismatch;
			 size_t minLenPrefix;
			 size_t minLenSuffix;
		    protected:

			virtual AbstractHit* createHit()
			    {
			    return new MismatchHit;
			    }
		    public:
			MisMatch(IndexedPositions* owner):PerfectMatch(owner),
			    maxMismatch(4),
			    minLenPrefix(5),
			    minLenSuffix(5)
			    {

			    }

		    virtual ~MisMatch() {}

		    virtual int compare2(IndexedPosition* pos,const CharSequence* seq,AbstractHit* hit)
			 {
			 FaidxPtr faidx= owner->genome->getFaidByTid(pos->chrom_id);
			 int64_t position=pos->position;
			 ChromosomeString chrom(owner->genome,faidx);
			 size_t seqidx=0;
			 int countMismatch=0;
			 //bool debug=strncasecmp( TEST_SEQ, seq,seqLen)==0;
			 //bool debug=true;
			 for(;;)
			      {
			      if(seqidx == seq->size())
				    {
				    //if(debug) cerr << "return 0\n";
				    return 0;
				    }
			       else if(position >= chrom.size())
				    {
				    return 1;
				    }

				int c1= normalize_base(chrom.at(position));
				int c2= normalize_base(seq->at(seqidx));
				//if(debug) cerr << "comparing["<< seqidx<<"/" << seq->size() <<"|"<< position <<"]"<<(char)c1 <<" / " << (char)c2 << endl;
				int i=c1-c2;
				if(i!=0)
				    {
				    if(seqidx < minLenPrefix ||
				       //seqidx+  minLenSuffix > seq->at(seqidx) ||
				       countMismatch== maxMismatch )
					{
					//if(debug) cerr << "break\n";
					return i;
					}
				    countMismatch++;
				    }
				seqidx++;
				position++;
				}
			 }

		    virtual int compare3(IndexedPosition* pos,const CharSequence* seq,AbstractHit* hit)
			 {
			 FaidxPtr faidx= owner->genome->getFaidByTid(pos->chrom_id);
			 ChromosomeString chrom(owner->genome,faidx);
			 size_t seqidx=seq->size()-1;
			 int64_t position=pos->position+seqidx;
			 int countMismatch=0;

			 for(;;)
				{
				int c1= normalize_base(chrom.at(position));
				int c2= normalize_base(seq->at(seqidx));
				//cerr << "comparing["<< seqidx<<"/" << seq->size() <<"|"<< position <<"]"<<(char)c1 <<" / " << (char)c2 << endl;
				int i=c1-c2;
				if(i!=0)
				    {
				    if(countMismatch== maxMismatch )
					{
					//cerr << "break\n";
					return i;
					}
				    countMismatch++;
				    }
				if(seqidx == 0)
				    {
				    //cerr << "return 0\n";
				    return 0;
				    }
			       else if(position==0)
				    {
				    return 1;
				    }
				seqidx--;
				position--;
				}
			 }


		    virtual int64_t find(const CharSequence* seq)
			{

			AbstractHit* hit=createHit();
			hit->seq=seq;
			hit->count=0;

			/* try left to right */
			if(seq->size()> minLenPrefix)
			    {
			    const SubStringSequence prefix(seq,0,this->minLenPrefix);
			    long n= lowerBound(&prefix);
			    while(n +prefix.size() <=  owner->size() )
				{
				owner->get(n,&(hit->position));
				if(compare(&(hit->position),&prefix)!=0) break;
				//cerr << "compare=0\n";
				if(this->compare2(&(hit->position),seq,hit)==0)
				    {
				    hit->count++;
				    hitFound(hit);
				    }
				++n;
				}
			    }
			/* try right to left */
			if(seq->size()> minLenSuffix)
			    {
			    const SubStringSequence suffix(seq,seq->size()-this->minLenSuffix,this->minLenSuffix);

			    long n= lowerBound(&suffix);
			    while(n +suffix.size() <=  owner->size() )
				{
				owner->get(n,&(hit->position));
				if(compare(&(hit->position),&suffix)!=0) break;
				hit->position.position -= seq->size()-suffix.size();
				//cerr << "compare=0\n";
				if(this->compare3(&(hit->position),seq,hit)==0)
				    {
				    //cerr << "OK rev"<< endl;
				    hit->count++;
				    hitFound(hit);
				    }
				++n;
				}
			    }
			int64_t count= hit->count;
			delete hit;
			return count;
			}

		    };




	    void dump()
		{
		IndexedPosition p;
		for(long n=0;n< size();++n)
		    {
		    get(n,&p);
		    FaidxPtr faidxptr=genome->getFaidByTid(p.chrom_id);
		    cout << (int)p.chrom_id <<":" << p.position << ":";
		    for(int i=0;i< 70 && p.position+i < faidxptr->len ;++i)
			{
			cout << genome->at(faidxptr,p.position+i);
			}

		    cout << "\n";
		    }
		}




    };



/*
static int default_hit_call_back_t(const IndexedPositions::HitDef* hitdef,void* userData)
    {
    FaidxPtr faidx=hitdef->indexed->getGenome()->getFaidByTid(hitdef->position.chrom_id);
    cout << "Match " << faidx->name << " at " << hitdef->position.position << endl;
    cout << " "; cout.write(hitdef->readSeq,hitdef->readLength); cout<< endl;
    cout << " ";
    for(size_t i=0;i< hitdef->readLength;++i)
	{
	cout << hitdef->indexed->getGenome()->at(faidx,hitdef->position.position+i);
	}
    cout << endl;
    return 0;
    }*/

#define UNSORTED_FILE "/tmp/jeter.idx"

static  int comparator(const IndexedPosition *a, const IndexedPosition *b,const void* userData)
    {
    IndexedGenome* genome=(IndexedGenome*)userData;


    const FaidxPtr faidx1=genome->getFaidByTid(a->chrom_id);
    const FaidxPtr faidx2=genome->getFaidByTid(b->chrom_id);

    int64_t n1=a->position;
    int64_t n2=b->position;

    for(;;)
	{
	if(n1 >= faidx1->len && n2 >= faidx2->len) return 0;
	else if(n2 >= faidx2->len) return 1;//YES
	else if(n1 >= faidx1->len) return -1;//YES ok for linux sort after dumping data
	int c1= normalize_base(genome->at(faidx1,n1));
	int c2= normalize_base(genome->at(faidx2,n2));

	int i=c1-c2;
	if(i!=0) return i;
	n1++;
	n2++;
	}

    }

static void indexGenome(const char* filename,const char *idxname)
    {
    IndexedGenome genome1(filename);
    FILE* out=safeFOpen(UNSORTED_FILE,"wb");
    int nRef=genome1.size();
    if(nRef>UCHAR_MAX)
	{
	cerr << "too many references in "<< filename << endl;
	exit(EXIT_FAILURE);
	}
    IndexedPosition data;
    for(int i=0;i< nRef;++i)
	{
	const FaidxPtr faidx=genome1.getFaidByTid(i);

	data.chrom_id=(unsigned char)i;
	for(int64_t n=0;n< faidx->len;++n)
	    {
	    if(normalize_base(genome1.at(faidx,n))=='N') continue;
	    data.position=n;
	    safeFWrite(&data,sizeof(IndexedPosition),1,out);
	    }
	}
    safeFFlush(out);
    fclose(out);

    /*
    cout << "AACTGCAACTCCAAAGCCACCCCTCACCCACTAGGATACCAACAAACCTACCCACCCTTAACAGTACATG\n";
    for(int i=0;i< 70;++i)
	{
	cout << genome1.at(genome1.getFaidByTid(0),16241+i);
	}
    cout << endl;*/

    MergeSort<IndexedPosition> merger(comparator);
    merger.sort(UNSORTED_FILE,idxname,&genome1);
    remove(UNSORTED_FILE);

    //IndexedPositions indexed(&genome1,idxname);
    //indexed.dump();
    }

static void findRead(const char* fasta,const char *idxname,const char* seq)
    {
    IndexedGenome genome(fasta);
    IndexedPositions indexed(&genome,idxname);
    IndexedPositions::PerfectMatch m(&indexed);
    CStringWapper str(seq);
    m.find(&str);
    }

struct FastQ
    {
    char* name;
    size_t len;
    char* seq;
    char* name2;
    char* qual;
    };



static void sampe(const char* fasta,const char *idxname,
	const char* fastQFile1,
	const char* fastQFile2
	)
    {
    IndexedGenome genome(fasta);
    IndexedPositions indexed(&genome,idxname);
    gzFile fastq1=safeGZOpen(fastQFile1,"r");
    gzFile fastq2=safeGZOpen(fastQFile2,"r");
    GZLineReader in1(fastq1);
    GZLineReader in2(fastq1);
    size_t nLine=0;
    for(;;)
	{
	nLine++;
	DefaultShortRead* reads[2];
	reads[0]=DefaultShortRead::read(&in1);
	reads[1]=DefaultShortRead::read(&in2);
	if(reads[0]==NULL || reads[1]==NULL) break;
	if(!(reads[0]!=NULL && reads[1]!=NULL))
	    {
	    throw std::runtime_error("not paired and files.");
	    }

	for(int i=0;i< 2;++i)
	    {
	    size_t len=reads[i]->size();
	    size_t offset=0;
	    while(len>0 && reads[i]->at(len-1)=='N')
		{
		len--;
		}
	    while(len>0 && offset< len &&  reads[i]->at(offset)=='N')
		{
		offset++;
		len--;
		}
	    if(len<10)
		{
		//cerr << "ignoring "<< reads[i].seq << endl;
		continue;
		}
	    TrimmedShortRead forward(reads[i],offset,len);
	    //CStringWapper forward(&trimmed,offset,len);
	    ReverseComplementString reverse(&forward);
	    bool found=false;
	    for(int strand=0;strand<2;++strand)
		{
		const CharSequence* shortReadSeq=(strand==0?(const CharSequence*)&forward:(const CharSequence*)&reverse);


		IndexedPositions::PerfectMatch m(&indexed);
		if(m.find(shortReadSeq)!=0)
		    {
		    found=true;
		    }


		if(!found)
		    {
		    IndexedPositions::MisMatch m2(&indexed);
		    if(m2.find(shortReadSeq)!=0)
			{
			found=true;
			}

		    }

		}
	    if(!found)
		{
		cerr << "cannot find:\n ";
		forward.print(cerr); cerr << "\n ";
		reverse.print(cerr); cerr << "\n";
		cerr << endl;

		}
	    else
		{
		cerr << "OK" << endl;
		}

	    }
	}

    gzclose(fastq1);
    gzclose(fastq2);
    }


int main(int argc,char **argv)
    {
    if(argc<2) return EXIT_FAILURE;
    if(argc==4 && strcmp(argv[1],"index")==0)
	{
	indexGenome(argv[2],argv[3]);
	}
    else if(argc==5 && strcmp(argv[1],"find")==0)
    	{
	findRead(argv[2],argv[3],argv[4]);
    	}
    else if(argc==6 && strcmp(argv[1],"sampe")==0)
	{
	sampe(argv[2],argv[3],argv[4],argv[5]);
	}
    else
	{
	cerr << "illegal arguments\n";
	}
    return 0;
    }
